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The spatial fluctuations of the extragalactic background light trace the total emission from 
all stars and galaxies in the Universe. A multi-wavelength study can be used to measure the 
integrated emission from first galaxies during reionization when the Universe was about 500 
million years old. Here we report arcminute-scale spatial fluctuations in one of the deepest 
sky surveys with the Hubble Space Telescope in five wavebands between 0.6 and 1.6 pm. We 
model-fit the angular power spectra of intensity fluctuation measurements to find the ultra¬ 
violet luminosity density of galaxies at z > 8 to be logpuv = 27. erg s“^ Hz~^ Mpc“® 
(Icr). This level of integrated light emission allows for a significant surface density of fainter 
primeval galaxies that are below the point source detection level in current surveys. 


Introduction 

The formation and early evolution of the first galaxies in the universe oeeurred some time after the 
dark ages, when the eoaleseenee of gravitationally bound masses formed in eomplex struetures, 
with a spatial distribution that ean be traeed baek to primordial overdensitie^iE. The ultraviolet 
(UV) photons from these first sourees initiated the reionization of the surrounding neutral medium, 
thus ending the dark ages and beginning the era of a transparent eosmos, which we are inereas- 
ingly familiar with today. The luminosity per unit volume of these UV photons at a rest wavelength 
around ISOOA (puv) during this reionization period is an important quantity to measure, as it traees 
the star formation and evolution of these ionizing sourees. The traditional method to measure the 
UV luminosity density of the universe, puv. during the epoeh of reionization, involves searehing 
for eandidate galaxies at z > 6 through their Lyman-dropout signature!^ and then eonstrueting 
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the luminosity function of those detected galaxies based on the observed number counts. This lu¬ 
minosity function is then extraploated to a fainter absolute magnitude and integrated in luminosity 
to calculate puv- 

There is a second way to quantify puv- This involves a measurement of the extragalactic 
background light (EBL) and, in particular, the angular power spectrum of the EBL intensity fluctu¬ 
ations. Because these intensity fluctuations are the result of emissions throughout the cosmic time, 
the signal we measure today is the sum of many different emission components, from nearby in our 
Galaxy to distant sources. If the integrated intensity from reionization can be reliably separated 
from that of foreground signals, we may be able to make an accounting of the total luminosity 
density of UV photons from reionization. Just as Eyman-dropout galaxies are detected in deep sky 
surveys, there is a way to achieve such a separation. Due to redshifting of the photons arising from 
sources present during reionization, their emission, as seen today, is expected to peak between 0.9 
and 1.1 /im. This assumes that the reionization occurred around 2 ; ~7 to 9, consistent with optical 
depth to electron scattering as measured by PlanclP. Due to absorption of ionizing UV photons, 
there is no contribution shortward of the redshifted Eyman break around 0.8 /irrl21121. Spatial fluc¬ 
tuations of the EBE centered around 1 /im thus provide the best mechanism to discriminate the 
signal generated by galaxies present during reionizatioiJiinil from those at lower redshifts, based 
on the strength of the drop-out signature in the fluctuations measured in different bands. 

There are existing measurements of the EBE flcutuations though their origin remain uncer¬ 
tain. This is mostly due to the fact that the previous measurements of EBE fluctuations have until 
now been limited to wavelengths greater than 1.1 /im, with the best measurements performed at 
3.6 /imP^. These studies have been interpreted with models involving populations of sources 
present during reionization at redshift z > 8, direct collapse and other primordial blackholes at 
z > 12 (Ref. [UnE)^ and with stellar emission from tidally stripped intergalactic stars residing in 
dark matter halos, or the “intra-halo light” (IHE^at 2 : < 3. The IHE is diffuse stars in dark matter 
halos due to galaxy mergers and tidal interactions. While the relative strengths of these various 
contributions are still unknown, we expect the signal from high-redshift galaxies to be separa¬ 
ble from low-redshift contributions, including those from faint nearby dwarf galaxie^, through a 
multi-wavelength fluctuation study spanning the 1 /rm range, including in the optical (A < 1 jum) 
and near-IR (A >1 jum) wavelengths. 

Here we present results from a multi-wavelength fluctuation study using data from the Hub¬ 
ble Space Telescope (HST) that span across the interesting wavelength range centered at 1 /im. 
Through models for multiple sources of intensity fluctuations, from diffuse Galactic Eight to pri¬ 
mordial faint galaxies, we are able to describe the five-band fluctuation measurements in optical 
and near-IR wavelengths to obtain a constraint on the UV luminosity density of galaxies present 
during reionization at a redshift above 8. We compare our measurement to existing constraints on 
the quantity and we also discuss implication of our measurement. 
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Results 

Fluctuation Power Spectra. We make use of imaging data from the Cosmic Assembly Near- 
Infrared Deep Extragalactic Legacy Surve 5 Ell^ (CANDELS), a legacy program of the Hubble 
Space Telescope (HST; see Eigure 1). Due to extensive data in the Hubble arch ive, we selected 
the southern area of the Great Observatories Origins Deep Survey (GOODS for the measure¬ 
ments (see the Methods section for details on our field selection process). This field contains HST 
observations with two instruments (Wide Eield Camera 3 and Advanced Camera for Surveys) that 
have the deepest and most continuous coverage. Our total dataset is comprised of observations that 
were taken between 2002 and 2012, with exposure times ranging between 180 seconds and 1469 
seconds per frame. We avoided the background gradients evident in the publicly available WEC3 
and ACS mosaics by creating our own self-calibrated mosaics using custom software!^, to produce 
120 square arcminute mosaics combined from 234 to 428 (depending on the passband) individ¬ 
ual flux-calibrated, flat-fielded (ELT) frames (Eigure 2). These mosaics are publicly available (see 
Supplementary Information Note 1). We mask stars and galaxies using an internally developed 
masking algorithm, facilitated by a public multi-wavelength catalog spanning from the ultraviolet 
to the mid-infrarecP^. The auto- (Eigure 3) and cross-power spectra (Eigure 4) are computed using 
standard Eourier Transform techniques on the masked images, which retain 53% of their pixels 
after masking. A number of corrections are performed on the power spectra. Details regarding 
these corrections can be found in the Methods section. 

Our measurements continue to show the significant excess in the fluctuation amplitude at 
30 arcsecond and larger angular scales, when compared to the clustering of faint, low-redshift 
galaxie^. The large-scale fluctuations correlate between filters (Eigure 4). The excess in the am¬ 
plitude of fluctuations relative to faint low-redshift galaxies is consistent with previous measure¬ 
ments at 3.6 Our HST-based power spectra probe deeper into the fluctuations and have 

shapes departing from the fluctuations measured with the GIBER sounding rocket experimenP^. 
Due to the shallow depth of the GIBER imaging data, the measured fluctuations there are domi¬ 
nated by the shot noise of the residual galaxies at the arcminute angular scales that we probe here 
with Hubble data (Eigure 5). 

At angular scales of tens of arcminutes and above, GIBER detected an up-turn in the fluctu¬ 
ations with an amplitude well above the level expected from instrumental systematics and residual 
flat-field errord^. However, as shown in Eigure 5, the combination of GIBER and Hubble fluctu¬ 
ations is consistent with a power-law clustering signal out to the largest angular scales probed by 
GIBER. If the power-law signal is of the form Ci oc the best-fit slope to combined GIBER and 
Hubble measurements at 1.6 fxm is a = —3.05 ± 0.07. This slope is consistent with Galactic dust, 
which in emission at 100 fim has a power-law of —2.89 ± 0.22 (Ref. At the largest angular 
scales we could be detecting interstellar light scattered off of Galactic dust or diffuse Galactic light 
(DGL). The overall amplitude of fluctuations we measure at 1.6 /rm is consistent with 10% DGL 
fluctuations at tens of arcminute angular scales, given the 100 fim surface brightness of GOODS-S 
of ~ 0.5 MJy/sr, and existing DGL intensity measurements (Eigure 6). Euture measurements in the 
optical wavelengths over a wider area are necessary to confirm the Galactic nature of fluctuations 
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at angular scales greater than one degree. 


Multi-Component Model. For our theoretieal interpretation, we invoke a model whieh involves 
four main eomponents: (a) intra-halo light (IHL) following Ref.l^, (b) diffuse galaetie light (DGL) 
due to interstellar dust-seattered light in our Galaxy, (e) low-redshift residual faint galaxie^, and 
(d) the high-redshift signal. We assume the flat ACDM model with VLm = 0.27, VLb = 0.046, 
as = 0.81, = 0.96 and h = 0.71 in our theoretieal modelin^^. We summarize the basie 

ingredients of our model now while providing referenees for further details. 

For IHL, we follow the model developed in Ref. The mean luminosity of the IHL at 
rest-frame wavelength A for a halo with mass M at 2 ; is deseribed as 

Lf^{M,z) = /ihl(M)L,^(M)(1 + ^)“F,(Ao/(1 + ;^)), (1) 

where Aq = A(1 -f 2 ;) is the observed wavelength, a is the power-law index whieh takes aeeount of 
the redshift-evolution effeet, and /ihl(^) is the IHL luminosity fraetion of the total halo, whieh 
takes the form 

/ M A^ 

/iHL(Af) = Aihl ] • (2) 

Here Ajhl is the amplitude faetor, Mq = 10^^ Mq, and f3 is the mass power index. In equation ([I]), 
L\^{M) = Lo{M)/Xp is the total halo luminosity at 2.2 /im and at 2 ; = 0, where Ap = 2.2 /im and 
Lq is given 

f M \ 

Lo(M)^5.64xlO-/r-(^^-^^^^ V (3) 

Here Hq = 70 /lyo kms“^Mpc“^ is the present Hubble eonstant. The Fx term is the IHL speetral 
energy distribution (SED), whieh ean transfer Lx^ to the other wavelengths and is normalized to 
be 1 at 2.2 /rm (see the diseussion in Ref. ^ for details). We assume the IHL SED to be the same 
as the SED of old elliptieal galaxies, whieh are eomprised of old red star^. 

The angular power speetrum of IHE fluetuations, G™'", is ealeulated via a halo model 
approaehf^ and involves both a 1-halo term assoeiated with spatial distribution inside a halo and a 
2-halo term involving elustering between halos. The details are provided in Ref. In the elus- 
tering ealeulation, we assume the IHE density profile follows the Navarro-Erenk-White (NEW) 
profileinil^. We set the maximum IHE redshift at z^ax = 6. The Mmin and Mmax are fixed to 
be 10® and 10^^ Mq/H, and the power-law index (3 is fixed to be 0.1 in this workpl. The IHE 
model is then deseribed with two parameters: Amh in equation ([^ and the power-law index a in 
equation Q. 

A simple test of IHE is to grow the souree mask and study how the fluetuation power spee¬ 
trum varies as a function of the mask size. However, we note that our model for IHE involves 
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clustering at large angular seales. That is, in our deseription IHL is not restrieted to regions near 
galaetie disks only. The dependenee of the power speetrum with mask size is studied in Refs.l^^Hl. 
These studies find that the fluetuations do not vary strongly with the mask radius, though sueh 
studies ignored the mode-eoupling effeets assoeiated with the mask as the mask radius is varied. 
As diseussed in Ref.^, in order to test IHL one has to grow the masking area to a faetor of 10 larger 
than the typieal mask radius used in the eurrent analyses of fluetuations. In the ease of Spitzer data, 
where we expeet fluetuations to be dominated by IHL, the relatively large 2 areseeond point spread 
funetion makes it elose to impossible to test IHL direetly with a varying mask. However, additional 
tests of IHL exist in the literature!^. These inelude eorrelations between artifieial halos and masked 
sourees, and eorrelations between masked sourees and foreground galaxies. Sueh tests have not 
ruled out the IHL eomponent. Furthermore, without sueh a eomponent we are not able to explain 
the fluetuations measured at wavelengths below 0.8 /im, as residual fiant galaxy elusterin^^is not 
adequate to explain the measurements. 

The DGL eomponent involves dust-seattered light and it is likely that the same dust is ob¬ 
served by IRAS at far-infrared wavelengths through thermal emission. The DGL eomponent was 
eonsidered in Ref.^ and an upper limit on the expeeted amplitude was ineluded based on the 
eross-eorrelation with 100 /im IRAS map of the same fields. The GIBER final model results fo- 
eussing on IHL to explain the fluetuations did not allow the DGL fluetuations amplitude to vary 
as a free parameter. With Hubble data, we find stronger evidenee for DGL or a DGL-like signal 
onee eombined with GIBER, and with an rms amplitude for fluetuations that is at least a faetor of 
3 larger than the upper limit used in Ref. ^ based on the eross-eorrelation with GIBER. Moreover, 
we find that the angular power speetrum is proportional to ~ over the degree seales measured 
from GIBER to tens of areminute seales of CANDEES measurements. So we model it with an 
amplitude faetor Adgl as 

= AoGLrl (4) 

To validate this DGE slope dependenee, we perform a linear fit in log spaee at low multipoles 
of the HST 1.6 /rm data simultaneously with the GIBER 1.6 /im data. We measure a slope of 
—3.05 ± 0.07, so the funetional form of our DGE model with is appropriate (Eigure 5). The 
power-law behavior of the DGE signal is eonsistent with Galaetie dust emission power speetra in 
far-infrared and sub-mm surveyP^, and dust polarization measurements in all-sky experiments like 
PlanekP^. We summarize a eomparison of our DGE intensity measurements to those of existing 
measurements as a funetion of wavelength in Eigure 5. 

The elustering of low-redshift faint galaxies at 2 ; < 5, where reliable luminosity funetions 
exist in the literature, is based on the detailed models in Ref. 1^. We follow the ealeulations pre¬ 
sented there to establish the expeeted level of low-redshift elustering, and the uneertainty of that 
expeetation at the depth to whieh we have masked foreground galaxies. Beeause the low-redshift 
luminosity funetions are not steep, unless there is a break or steepening in the luminosity funetion 
— whieh is not supported by the halo model — these low-redshift populations do not dominate the 
elustering we have measured. Given that our fluetuation measurements reaeh the deepest depths 
provided by both Hubble and Spitzer/IRAC, it is also unlikely that populations sueh as extreme red 
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galaxies at z < 2 are responsible for the measured fluetuations. If there are populations at low red- 
shift responsible for the SED of the fluetuations we measure, they would need to have individual 
SEDs that are eonsistent with a sharp break redshifted between 0.8 and 1.25 /rm. While fluetuation 
measurements in just two bands eannot separate galaxies that have redshifted 4000 A break, or 
galaxies that have redshifted Eyman-a break between those two bands, with five bands we have 
adequate knowledge on the SED of fluetuations, and the shape of the elustering over two deeades 
in angular seales to separate high- 2 ; galaxies from low- 2 ; faint interlopers. The low- 2 ; interlopers 
are also likely eaptured by our IHE model as we eannot distinguish between diffuse stars and faint, 
dwarf galaxies that happen to be a satellite of a large dark matter halo in our modeling deseription. 


Galaxies During Reionization. The final and eritieal eomponent in our model is the signal from 
2 ; > 6. We break this signal into two redshift intervals given the plaeement of the five ACS and 
WEC3 bands, based on the Eyman-dropout signal that moves aeross these bands. In partieular, we 
eonsider 8 < 2 ; < 13 and 6 < 2 ; < 8 as the two windows. As we diseuss later, given the availability 
of SER density measurements in the 6 < < 8 interval, we mostly allow the signal in that redshift 

interval to be eonstrained by the existing data, and model-fit independently the SER density in the 
higher-redshift interval. We do not have a strong independent eonstraint on the 6 < .^ < 8 signal 
since it is only a Eyman-dropout in our shortest-wavelength band at 0.6 /im. This allows a better 
separation of the 8 < 2 ; < 13 from the rest of the signals discussed above. To measure 6 < 2 ; < 8 
independently, we would need at least one more band below 0.6 /im. The signal from 8 < 2 ; < 15 
disappears from the three optical bands (0.6 to 0.85 /xm) and is present in the two IR bands at 1.25 
and 1.6 /xm. 


To model the high-redshift signals, we adopt an analytic modef^I^ based on the work of 
Ref. 1^. It involves a combination of two separate classifications of stars — moderate-metallicity, 
second-generation or later stars (PopII), and the first generation of stars ever formed in the Uni¬ 
verse, hence zero metallicity (PopIII). These are modeled with Salpetei^^and Earsorpl initial mass 
functions (IMEs) for PopII and PopIII stars, respectively. The calculation related to direct stel¬ 
lar emission and the associated nebular lines, including especially Eyman-a emission, follows the 
work of Eemandez & KomatsiP^. The total integrated intensity from Zj^in < z < z^i^ is 





C iy{z)ju{z) 

^H{z) {1 + zy 


(5) 


where u{z) = {1 + z)^. The comoving specific emissivity, as a function of the frequency is 
composed of both PopII and PopIII stars with an assumed z-dependent fraction as discussed in 
Ref. with the form given by 


fp{z) 


1 -f erf 


10 


(Jr) 


( 6 ) 


with (Tp = 0.5. The model thus assumes most of the halos have PopIII stars at z > 10 while PopII 
stars dominate at redshifts lower than that. 
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There are a number of theoretieal parameters related to this model, espeeially the eseape 
fraetion of the Lyman-a photons /esc, the star-formation effieieney denoting the fraetion of the 
baryons eonverted to stars in high-redshift dark matter halos, or /*, and the minimum halo mass to 
host galaxies, or Mmin. The overall quality of the data is sueh that we are not able to independently 
eonstrain all of the parameters related to the high-redshift intensity fiuetuation signal. Moreover 
these parameters are degenerate with eaeh other (i.e., ehanging /* ean be eompensated by a ehange 
in Mmin for example). Thus given that we do not have the ability to eonstrain multiple parameters, 
we simply model-fit a single parameter Ahigh-z that seales the overall amplitude from the default 
model, interpret that sealing through a variation in /*, and subsequently eonvert that to a eonstraint 
on the SFRD. We fix our default model to a basie set of parameters, and assume /esc = 0.2, 
Mmin = 5 X 10^ Mq, and /* = 0.03. The resulting optieal depth to reionization of this default 
model is 0.07, eonsistent with the optieal depth measured by PlanelP. Among all these parameters, 
the most signifieant ehange (over the angular seales on whieh we measure the fluetuations) eomes 
effeetively from /*, or the overall normalization of jy{z), given that it is direetly proportional to 
/*. This ean in turn be translated to a direet eonstraint on the SFRD, '/(^), sinee with /* we are 
measuring the integral of the halo mass funetion sueh that 


where dn/dM is the halo mass funetiorl^. 




nT) 

dMM—(M, z) 
dM^ ’ ^ 


(7) 


Finally, to ealeulate the angular power speetrum of fluetuations, we also need to assign galax¬ 
ies and satellites to dark matter halos. For that we make use of the halo modeP^. We make use 
of the same oeeupation number distribution as in Ref. ^ where the eentral and satellite galaxies 
are defined following Ref.^. However, departing from the low-redshift galaxy models, we take a 
steep slope for the satellite eounts in galaxies with = 1-5. The low-redshift galaxy elustering 
and luminosity funetions are eonsistent with ~ 1 (Ref. 1^, but sueh a value does not reproduee 
the steep faint-end slopes of the LEG luminosity funetionP^. Sueh a high slope for the satellites 
also boost the non-linear elustering or the 1-halo term of the fluetuations. We do not have the 
ability to independently eonstrain the slope of satellites from our fiuetuation measurements. In 
the future a joint analysis of fluetuations and LEG luminosity funetions may provide additional 
information on the parameters of the galaxy distribution that is responsible for fluetuations. It may 
also be that the models ean be improved with additional external information, sueh as the optieal 
depth to reionization. We also note that other sourees at high redshift inelude direet eollapse blaek 
holes (DCEH^i^, but we do not explieitly aeeount for them here as the existing DCEH model is 
finely tuned to mateh Spitzer fluetuations, and the low signal-to-noise ratio of the Chandra-Spitzer 
eross-eorrelation results in them residing primarily at z > 12. DCEHs at sueh high redshifts will 
not eontribute to Hubble fluetuations. 


Finally, at smaller angular seales, the shot noise dominates the optieal and IR baekground 
intensity fiuetuation. Sinee it is seale-independent, we set it as a free variable parameterized as 

C'f = Ashot. (8) 
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This noise term in the fluctuation power spectrum arises because of the Poisson behavior of the 
galaxies at small angular scales, a product of the finite number of galaxies. Our measured shot 
noise comes from a combination of the unmasked, faint low-redshift dwarf galaxies, and the high- 
redshift population. We do not use the information related to the shot noise in our models but 
instead treat it as a free independent parameter, since we cannot separate the high-redshift shot 
noise from the shot noise produced by faint, low-redshift dwarf galaxies. Here we focus mainly on 
the clustering at tens of arcseconds and larger angular scales to constrain SFRD during reionization. 
In the future, with either a precise model for the low-redshift galaxies or a model for high-redshift 
galaxies that determines their expected number counts as a function of the free parameters such 
as Mmin and /*, it may be possible to separate the overall shot noise associated with reionization 
sources from that of the low-redshift faint galaxies. If this is the case then it might also be possible 
to improve the overall constraints on the high-redshift population. It may also be that, under an 
improved model, shot noise may end up providing complementary information to galaxy clustering 
to break certain degeneracies in model parameters. 

Our overall model for the optical and infrared background fluctuations is 

( + Cf -f F125W and above 

Ce = } C-IHL + C'^GL ^low-z ^shot ^ q6<z<8 F775W and F850LP (9) 

[ F606W 

Given that we are not able to constrain the the amplitude of given the degeneracies with 

the parameters involving the IHL model, and the fact that we only have a single band below it, we 
set based on the default prediction of our model, but allow the overall amplitude 

to vary such that it uniformly samples the SFRD between [0.003,0.2] Mq yr“^ Mpc“^. The range 
is fully consistent with the existing measurements on the SFRD between z = 6 and 8 (Ref. 1^®). 
Our constraint on is mostly independent of this parameter since we can safely constrain 

the Lyman-dropout signal between 0.8 and 1.25 /rm with our existing data. 

We also included the CIBER^data at 1.1 and 1.6 /rm and Spitzei^i^ data at 3.6 /im in our 
fitting process. When compared to the Hubble data at 1.25 and 1.6 /im, we find the GIBER data are 
likely dominated by the emission from a DGE-like signal at large angular scales, and low-^ faint 
galaxies at 2 ; < 5 at small angular scales (Eigure 5). Eor the fluctuations from faint, low-;^ galaxies, 
we adopt a model of residual galaxies which is derived from the observations of the luminosity 
function for different near-IR bandd^. This model already includes the shot noise term, and we 
add a scale factor fiow-z to vary the low- 2 ; angular power spectrum, in 1 a uncertainty. Eor 

the DGE component, we use the of Hubble data at 1.25 and 1.6 /xm to fit the GIBER data at 
l.I and 1.6 /im. 

We perform joint fits for Hubble, GIBER and Spitzer data with the Markov Chain Monte 
Carlo (MCMC) method. The Metropolis-Hastings algorithm is used to find the probability of ac¬ 
ceptance of a new MCMC chain poin^i^EIl. We estimate the likelihood function as G cc exp(— x^/2). 
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where is given by 


X 


2 


^^obs _ ^th^2 

;;2 


( 10 ) 


Here and Cf^ are the observed and theoretieal angular power speetra for HST, Spitzer or 
CIBER data, respeetively. at is the error for eaeh data point at I, and is the number of data 
points. The total of HST, Spitzer and CIBER is Xtot = xIst + Xciber + Xspitzer- 


We assume a flat prior probability distribution for the free parameters; see Table 1 for prior 
information. Both Adgl, vary as independent parameters for eaeh band. Both the Adgl 
and Cf'°^ parameters are six-fold, with one for eaeh HST band and for Spitzer/IRAC 3.6 fim. (we 
eombined the two CIBER bands with two of the HST bands). We have two parameters for IHE 
and one parameter for the normalization of the reionization galaxies with A 8 < 2 <i 3 - We have two 
more parameters that we vary, A 6 < 2<8 and /low-z- We set a uniform prior on Aq^z<8 in the SERB 
following the existing measurements to be between 0.003 and 0.2 M© yr“^ Mpe“^. We also set a 
uniform prior on /low-z over a reasonable range of models to aeeount for the overall uneertainty in 
the models of Ref. ^ to deseribe the z < 5 faint galaxy elustering at the same masking depth as our 
measurements. We marginalize over both and fiow-z as well as all other parameters when 

quoting results for an individual parameter. We have a total of 14 free parameters in our MCMC 
fitting proeedure that we extraet from the data. Among these parameters, 12 of them simply de¬ 
seribe the small and large angular seale fluetuations in eaeh of the bands we have performed the 
measurements. These parameters are summarized in Table 1. We generate twenty MCMC ehains, 
where eaeh ehain eontains about 100,000 points after eonvergenee. After thinning the ehains, we 
merge all ehains and eolleet about 10,000 points for illustrating the probability distributions of the 
parameters. Contour maps for eaeh of the fitted model parameters are shown in Eigure 8 . Our best- 
fit model with 14 free parameters have a minimum value of 278 for a total degrees of freedom 
of A^dof = 104. 

Discussion 

Our results are summarized in Eigure 3, where we show the best-fit model eurves. While the 
dominant eontribution to the exeess fluetuations eomes from DGE at £ < 10^, at intermediate 
seales we find the IHE and reionization eontributions to be roughly eomparable. In Eigure 6 we 
show the rms fluetuation amplitude at ~ 5 areminute angular seales over the interval 10000 < i < 
30000. We find a speetral energy distribution that is eonsistent with Rayleigh-Jeans (RJ) from 4.5 to 
2.4 yum, but diverges between 2.4 and 1.6 fxm, and even more rapidly between 1.25 and 0.85 /im. 
The fluetuations ean be explained with a eombination of IHE and high-redshift galaxies. The 
residual low-z galaxy signal is small but non-negligible. We find that it is mostly degenerate with 
IHE, espeeially if we allow its amplitude to vary more freely than the range allowed by the existing 
models based on 2 ; < 5 galaxy luminosity funetion^. Thus modeling uneertainities related to 
the I 0 W- 2 ; galaxy eonfusion do not eontaminate our statements about reionization. Assuming the 
existing low- 2 ; galaxy modeP^, the best-fit model is sueh that the IHE intensity peaks at lower 
redshifts with deereasing wavelength (Eigure 7). At 3.6 /rm, the IHE signal is assoeiated with 
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galaxies at z ~ 1, while at 0.6 /xm over 80% of the signal is associated with galaxies at z < 0.5. 
The total intensities are 0.13]'^o;o5, 0.23^0;}^, 0.27loj3, 0.45lo;24, and O.hd^gJ^ nW m“^ sr“^ at 
0.60, 0.77, 0.85, 1.25 and 1.6 /xm, respectively. We find that the implied IHL intensities at 1.25 and 
1.6 /xm are a factor of 10 lower than the implied IHL intensities for a model of GIBER fluctuations 
with IHL alone. The difference is due to the GIBER model that only included IHL and ignored the 
presence of DGL. 

The drop in the fluctuation amplitude from 1.25/1.6 /xm to 0.85 /xm allows for a signal from 
reionization, but the presence of fluctuations at shorter wavelengths, such as 0.6 /xm, rules out a 
scenario in which reionization sources are the sole explanation for the fluctuations at wavelengths 
at 1 /xm and above. The 3.6 /xm and X-ray cross-correlatioiJ^ was explained with primordial direct 
collapse blackholes at 2 ; > 12 (Ref.^^. In our multi-component model we are able to account for 
the presence of fluctuations at short wavelengths with IHL, DGL and faint low-redshift galaxies, 
while a combination of those components and high-redshift galaxies are preferred to account for 
fluctuations at 1.25 and 1.6 /xm. The high- 2 ; signal is modeled following the calculations in Ref.^. 
The signal has an overall amplitude scaling that is related to the star-formation rate during reion¬ 
ization. The bright end of the counts are normalized to existing luminosity function measurements, 
and the faint-end of the luminosity functions to have a steeper slope than measured with counts ex¬ 
tending down to arbitrarily low luminosities. In order to test whether a component at high redshift 
is required to explain the measurements, we also re-ran the MGMG model fits but with Ahigh-z 
fixed at 0. In this case our best-fit model with 13 free parameters has a minimum value of 283 
for a total degrees of freedom of iVdof = 105. The difference in the best-fit values with and 
without a model for high-redshift galaxies suggests a p-value of 0.025. This is consistent with the 
2cr to 3cr detection of 8 < 2 ; < 13 signal in the fluctuations (Ligure 6). 

With multi-wavelength measurements extending down to the optical, we are now able to 
constrain the amplitude of that signal with a model that also accounts for low-redshift sources in 
a consistent manner. This improves over previous qualitative arguments that have been made, or 
models involving high-redshift sources alone that have been presented, for the presence of a signal 
from reionization in the IR background fiuctuation^i^lllllll. in our models, the total intensity arising 
from all galaxies at 2 ; > 6 is log = —0.32 ±0.12 in units of nW m“^ sr“^ at 1.6 /xm. At 1.6 /xm 
the intensity from high-redshift sources is dominated by 2 ; > 8 galaxies, while at 0.85 /xm we find 
an intensity logz/J^, = —0.75 ± 0.05 in units of nW m“^ sr“^ for 6 < 2 ; < 8 galaxies. The total 
intensity from 2 : > 8 galaxies in the 1.6 /xm band is comparable to the IHL intensity at the same 
wavelength (Ligure 7). However, at 3.6 /xm, the IHL signal is a factor of about 5 times brighter 
than the z > 8 galaxies. At 1.6 /xm the total of the IHL, high-z, and integrated galaxy lighP^of 
IO.OI 4 8 m“^ sr“^is comparable to the EBL intensity inferred by gamma-ray absorption datcP^ 

of 15 ± 2(stat) ± 3(sys). 

Using the best-fit model and uncertainties as determined by MGMG model fits, we also 
convert the A 8 <^<i 3 constraint to a measure of the luminosity density of the universe at 2 : > 8 
(Ligure 9). The resulting constraint is log puv = 27 in units of erg s~^ Hz“^ Mpc“^ at (1 a). 
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As shown in Figure 9, the 68% confidence level constraint on puv is higher than the existing results 
from Lyman drop-out galaxy surveys during reionization at z > 8 (Refs.15313), and especially at 
10 (Ref. SI). At the 95% confidence level, our measurement is fully consistent with the existing 
results at 2 ; ~ 10 (Ref. SSI). Our constraint allows for the possibility that a substantial fraction 
of the UV photons from the reionization era is coming from fainter sources at depths well below 
the detection threshold of existing Lyman dropout surveys, as is indeed anticipated from the steep 
measured slopes of the UV luminosity functions from detected galaxies. Despite their lack of 
detections in the deepest surveys with HST, the majority of the faint sources responsible for both 
fluctuations and reionization should be detectable in deep surveys with JWST centered at 1 fxm. 

Methods 

Field selection. In order to obtain angular power spectra over large angular scales, individual ex¬ 
posures must be combined into one or more mosaiced images. We generate our own mosaics using 
the self-calibration technique!^ (SelfCal), instead of using the publicly available mosaiced images 
produced by astrodrizzle (available at http://candels.ucolick.org/data_access/Latest_Release.html). 
Foreground emissions, predominately that of Zodiacal light, are particularly pernicious at infrared 
wavelengthd^, so care must be taken when producing mosaics which combine observations taken 
at different times, especially at the WFC3/IR channels. Offsets between frames will lead to a fic¬ 
titious anisotropy signal if care is not taken to properly model and remove those offsets, which is 
what SelfCal was designed to do. 

Although the CANDELS observations cover multiple fields, only the two deep fields—the 
Great Observatories Origins Deep Survey-South and -North (GOODS-S and GOODS-Nj^^have 
sufficient overlap between frames to perform a self-calibration. The GOODS-N dataset has a larger 
number of frames with clear overall offsets resulting from scattered light than does GOODS-S. 
Therefore we have restricted our analysis to GOODS-S, which has an area of approximately 120 
square arcminutes. The wider CANDELS fields are composed of much poorer tile patterns, as 
can be seen in Fig. 18 of Ref. which are a significant draw-back for fluctuation studies, since 
one cannot calibrate the full mosaic to a consistent background level without introducing artificial 
gradients to the background intensity. Such dithering patterns were pursued by CANDELS to 
maximize the total area covered with WEC3. The increase in area is of benefit to studies that aim 
to detect rare galaxies, such as Lyman-break galaxy (LBGs) at z > 5. 

Initial data reduction and map-making. In addition to the data collected by CANDELS, the 
GOODS-S field has a wealth of HST archive data, publicly available on the Barbara A. Mikul- 
ski Archive for Space Telescopes (MAST; located at https://archive.stsci.edu/hst/search.php). We 
assembled our own collection of calibrated, fiat-fielded (ELT) frames from the MAST archive, 
comprised of some or all of the data from ten different HST proposalsI^iH^^ini. These data are 
also supplemented by the Early Release Science observationd^(ERS). The tile patterns of these 
observations can be found in Eigure 1. 
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In addition to selecting frames with a favorable tile pattern appropriate for self-calibration, 
we also had to take two additional potential issues into account. After the replacement of the ACS 
CCD Electronics Box during the fourth Hubble servicing mission (SM4), ACS imaging data are 
plagued with horizontal striping dominated by 1/f noise. Furthermore, ACS frames have a tendency 
to introduce a Moire pattern (correlated noise) when the pixel scale is modified in a low signal-to- 
noise area. Both of these characteristics can potentially contaminate an angular power spectrum 
to such an extent that the systematics dominate the measurement. With simulations, we found that 
with an increased number of ACS frames taken with varying position angles effectively removes 
the Moire pattern upon repixelization, and the bias-striping issue is ameliorated by simply omitting 
a large percentage of post-SM4 frames. Any given collection of ACS frames we used contained 
< 27% post-SM4 frames. 

Our MAST archive data were initially reduced with PyRAF version 2.1.1. MAST queries are 
reprocessed “on-the-fiy”, which entails using the most recent calibration files. Thus the FFT frames 
we retrieved from the archive had standard calibrations of bias and dark frame subtraction, along 
with flat-field correction, already performed. We identified cosmic rays in the the FFT frames with 
the CRCLEAN PyRAF module; sub-arcsecond astrometric alignment against the publicly available 
CANDEFS mosaic^was achieved with tweakreg. All ACS data were charge transfer efficiency 
(CTE) corrected, and post-SM4 ACS frames were de-striped prior to CTE correction. 

We generate mosaics from the reduced FFT frames using the same SelfCal model as in Ref.^^, 
an example where HST data has already been self-calibrated; details of the model can be found 
there. We de-weight bad pixels and cosmic rays, and iterate three times in order to find a SelfCal 
solution. Our input FFT frames are geometrically distorted with a pixel size of 0."0498 x 0."0502 
forACS,andO."1354 X 0."1210 for WFC3. We remove the distortion in the map making procedure 
and produce mosaics with a slightly larger pixel size of 0."140 (geometrically square). 

Note that the algorithm we use is same as the one used to generate self-calibrated maps of 
IRAC in the Spitzer fluctuation studie^f^HHH. The same method was implemented by the Her- 
schel SPIRE Instrument Science team to generate wide area mosaics of the Herschel-SPIRE data, 
resulting in far-infrared fluctuation^. The algorithm originates from the time of FIRAI^ and has 
wide applications. In the future we expect it will be used to combine frames and produce stable 
wide-area mosaics from JWST, Euclid, and WFIRST, among others. 

We generate two maps per band so we can use the differences and sums to study systematics 
and noise biases, as was done in previous studied^. The data are sorted by observation date and 
every other FFT frame was used for each half map. HST data generally have two or more exposures 
per pointing, so this results in two maps per band of the same or similar dither pattern and exposure 
time per pointing. One map from each band can be found in Figure 2. Multiple maps of the same 
band enable us to do cross-correlations, which ensures a removal of uncorrelated noise in the auto- 

'http://candels.ucolick.org/data_access/GOODS-S.html 
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spectra. This jack-knife process is similar to all other analyses related to large-scale structure and 
CMB angular power spectra from maps. 

Generation of resolved source mask. We utilized existing multi-wavelength catalogs of de¬ 
tected sources from the ultraviolet to mid-infrared (CITO/MOSAIC, VLT/VIMOS, VLT/ISAAC, 
VLT/HAWK-I, and /IRACj^. In addition, all the sources from the CANDELS, HUDF and 

ERS surveys (E435W, E606W, E775W, E814W, E850EP, E098M, E105W, E125W, and E160W) 
are also present in the catalog. The 50% completeness limit for F160W in the catalog is ttiab = 
25.9, 26.6 and 28.1 for the CANDEES wide, deep and HUDE regions respectively; the 5 a limiting 
magnitudes are 27.4, 28.2 and 29.1. Eor each source detected by any of the aforementioned in¬ 
struments, we apply an elliptical mask with parameters corresponding to the S EXTRACTOR Kron 
elliptical aperture. This catalog simplifies our masking procedure and ensures we are masking 
sources detected at other wavelengths, which may otherwise be undetected in our five bands. 

In addition to the source mask generated from sources detected by other instruments, we 
also generate our own internal masks. We run SExtractor on a coadded map in each of our 
five filters and apply the same elliptical masking procedure to incorporate the shapes of sources. 
Next, we take the union of all five internal source masks, plus the source mask we made from 
the pre-existing catalog. After applying this union mask to each band, we clip 5 a outliers and 
visually inspect each masked map. Any residual sources are masked by hand. We verified that all 
sources detected above Sa in any of the bands, including deep IRAC data at 3.6 /rm, are masked. 
This process yields 53% of the pixels unmasked for the EFT computation. Note that tests can be 
performed which expand and shrink the source mask to further test the IHE model (Ref 1^. 

Absolute flux calibration. SelfCal achieves relative calibration between frames, so in general, the 
absolute flux calibration, or gain, of SelfCal output maps needs to be determined from a standard 
flux reference. The multi-wavelength catalog used in the masking procedure is in principle a good 
enough reference, however the photometry in the public CANDEES source catalogs was obtained 
with a private version of S EXTRACTOR, and has aperture corrections applied to the flux densities 
which were extracted from PSE-matched maps. We instead generate internal catalogs from public 
CANDEES MultiDrizzled mosaics using the same procedure as we used on our mosaics. 

In each band, we re-pixelize the MultiDrizzle maps to our SelfCal pixel scale, and perform 
source extraction with SExtractor on both our mosaics and the MultiDrizzle mosaics. We use 
the same parameter files for all the source extractions. We then astrometric ally match the resultant 
catalogs in each band and keep those sources that are common within a radius of 0".l (our ELT 
frames were aligned in TWEAKREG with MultiDrizzle mosaics as the astrometric reference, so 
our astrometry is similar to the MultiDrizzle maps at sub-arcsecond scales). Counts in e—/s are 
converted to /iJy using the current HST magnitude zero pointP^. The calibration introduces a 
4 — 5% error which is propagated into our final error bars. 
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Power spectrum evaluation. All the statistical information contained in any one of our maps is 
summarized by its angular power spectrum, Ci, which is just the variance of the a^m’s. We used 
standard Fast Fourier Transform (FFT) techniques to estimate the C^’s of our masked map^i^E^E^I^. 
The mosaics can be seen in Figure 2, both in real space and Fourier space. For each band we have 
two half maps, A and B. To measure the inherent noise in the data, we compute the noise power 
spectrum as the auto power spectrum of (A—B)/2. To measure the raw auto-spectrum, we compute 
the cross spectrum of the two half maps, AxB, which eliminates any uncorrelated noise in our 
power spectrum estimate. In general, the standard deviation at each multipole is 


5C. 


/sky(2£+l)A£ 


//oauto I /onoise\ 

-h J, 


( 11 ) 


where A£ is the bin width for the given Ci, and fgky is the fractional sky coverage from all the 
pixels used in the FFT (excluding zeros). In our case, we have some error associated with the 
absolute calibration of the maps, so we take the total error budget of our raw power spectra as 
the quadratic sum of the calibration errors with the variance in equation ( [TT] ). The final measured 
auto-spectra and associated errors can be found in Supplementary Table 1. 


To account for the effects that the source mask, tile pattern and finite beam size introduce 
into the power spectrum, we employ the correction techniques of the MASTER algorithmP^, and 
closely follow the implementation procedures explained in Section 4, 5 and 6 of the Supplementary 
Information (SI) of Ref. ^ (Including SI Figure 1 of Ref.^. Among the procedures listed there, 
we have only slightly modified the way we generate our transfer function, T(£). In addition to 
adding instrumental noise (step 2 in Section 6 of the SI of Ref.^^, we also add an offset to each 
tile equal to the median of the given FLT frame. This additional step should in principle be a good 
indicator of how well SelfCal is performing in offset removal, unique to each observation. Transfer 
functions for each of our five bands can be found in Supplementary Figure 2; measurements of the 
beam transfer function can also be found in Supplementary Figure 2. 
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Figure 1 | Summary of tile patterns and their data archive identifications, a, proposal ID’s 
for each filter in the GOODS-S field. The ACS and WFC3 rows show the proposals which are 
common between all the bands in each instrument. For each proposal we did not necessarily use 
all the frames, specifically those from deep surveys. Also show are the tiling patterns for all the 
bands: 0.606 /xm (F606W; b), 0.775 /xm (F775W; c), 0.850 /xm (F850LP; d), 1.25 /xm (F125W; 
e) and 1.60 /xm (F160W; f). The units of the tile pattern figures are Logio(N+l), where N is the 
number of frames overlapping. The dashed white line indicates the cropped region where the 
fluctuation analysis was performed. 
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Figure 2 | Self-Calibrated Mosaics, a, GOODS-S SelfCal mosaics for each band. We astromet- 
rieally align eaeh map and erop the outer regions so as to inelude only seetions that have been 
observed in all five bands. The units of the maps are nW m~^ sr“^. b, the same as a exeept 
with the souree mask applied, c shows the fast fourier transform (FFT) of eaeh of the maps in b, 
whieh is what is used to measure the angular power speetrum. This is plotted in Fourier spaee as 
a funetion of modes and iy. The FFT of eaeh map is struetureless and eontains only Gaussian 
noise, whieh is indieative of high quality mosaies. By definition, the units of the FFT are the same 
as the units of the map. Eaeh eolumn is filter speeifie, plotted as 0.606 /im, 0.775 /im, 0.850 /im, 
1.25 /im and 1.60 /rm. 
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Figure 3 | Angular power spectra of optical and near-infrared background intensity fluctu¬ 
ations. a: Multi-wavelength auto power spectra of optical to near-IR intensity fluctuations in the 
GOODS-South field using Hubble Space Telescope data (see Online Methods for data selection 
details). The error bars are calculated by adding in quadrature the errors from the beam transfer 
function, map-making transfer function and calibration errors, to the standard deviation at each 
multipole, bCi, described in Equation 11. Thus the 1 a uncertainties account for all sources of 
noise and error, including map-making, calibration, detector noise, and cosmic variance associated 
with finite size of the survey. We show the best-fit model which makes use of four components: 
(a) z > % high-redshift galaxies, (b) intra-halo light (IHL)[i^, (c) faint low-redshift galaxied^, and 
(d) diffuse Galactic light. At 1.25 and 1.6 /im, the best-fit high-redshift galaxy signal is shown 
as dashed lines. The signal is zero in the optical bands. We show the upper limit (denoted by 
a downard facing arrow) of fluctuations generated by 6 < 2 ; < 8 galaxies as a dot-dashed line. 
Fluctuation power spectra and the best-fit models with la error bounds for the model components 
are shown at 0.775 /im in b and 1.60 /xm in c. The dominant model contributors to the total 
power spectrum are DGL at low multipoles, or angular scales greater than a few arcminutes, IHL 
at intermediate multipoles corresponding to angular scales of about an arcminute, and shot noise 
associated with faint low-redshift dwarf galaxies dominating the high multipoles or sub-arcminute 
angular scales. The clustering signal of low-z galaxies is more than an order of magnitude below 
the lower limit plotted here, thus we did not include a low-;^ component in our modeling. 
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Figure 4 | Angular cross- power spectra of optical and near-infrared background. Ten 

cross-correlations between the HST bands. Excess signal is detected in the cross-correlations. 
The error bars are 1 a uncertainties which are calculated in a similar way as in Figure 3, which 
accounts for all sources of noise and error, including map-making, calibration, detector noise, 
and cosmic variance. However, the noise power spectra for the cross-correlations are calculated 
slightly differently. For each filter we have two maps, so for each cross-correlation between bands 
we have four maps (label them A and B for the first filter, and C and D for the second). This 
enables us to generate a noise power spectrum by computing (A-B) x (C-D), as opposed to taking 
the auto-spectrum of (A-B) for the auto-correlations. The first row corresponds to all correlations 
with the 0.606 /rm band, the second for all correlations with the 0.775 /rm band not found in the 
first row, the third row corresponds to all correlations with the 0.850 /rm band not found in any 
of the preceding rows, and the last row corresponds to correlations at 1.25 /rm. The columns 
similarly increase in wavelength as you move across the page. 
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Figure 5 | Various auto-spectra and the spectral energy distribution of diffuse galactic light. 
Left, CANDELS corrected power spectra plotted against the CIBERp^, Spitzei^ and NICMOSp^ 
measurements (see also Ref. ® for more recent NICMOS measurements). The power spectrum 
resulting from the NICMOS analysis was measured from a MultiDrizzle map and has not been 
corrected for the transfer function and mode-coupling matrix resulting from source masking as 
discussed in our Methods section. Therefore we show it as a comparison but do not use it in our 
modeling. The error bars are 1 a uncertainties that account for all sources of noise and error, 
including map-making, calibration, detector noise, and cosmic variance associated with finite size 
of the survey. Right: Optical and infrared diffuse galactic light (DGE) spectrum. The CANDEES 
points are taken from the DGE model components at 10^ < i < 3 X 10^, and the GIBER points 
are taken directly from Eig 2. of Ref.® where they subtract off the shot noise component from their 
data. The galactic latitude for the optical points are |6| ~ 39°, 32°, 41°, 40° for the points labeled 
WitP^, Pale)^^, lenakapl and Guhathakurtcpl. The BrandP^ points are modeled over the full sky. 
GOODS-S is at a galactic latitude of |6| = 54°. 
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Figure 6 | Spectral energy distribution of optical and infrared fluctuations at arcminute an¬ 
gular scale. The Hubble/CANDELS points are averaged over 10^ < i < 3 x 10"^, with the 
best-fit shot noise and DGL eomponents subtraeted. Our model fits for the high-redshift and IHL 
eomponents, with their 1 a bounds, are shown as the filled regions. The errors here are propagated 
from the errors on the auto speetrum at the same i range. The light blue region shows the 1 a 
eonfidenee bound for the IHL eomponent when we use only the HST data in our model fitting; the 
dark blue region shows the 1 a eonfidenee bound for the IHL eomponent when use both the HST 
and Spitzer IRAC data in our model fitting. The light red eolored region signifies the 1 a error 
bound for the high-redshift model eomponent. The dashed line corresponds to the 1 a bound for 
the low-redshift component. The Spitzei^ and AKARp^ data are taken from previous measure¬ 
ments at £ = 3000. Note the spectral dependence difference between the high-redshift signal and 
IHL. Below 0.8 /im we do not expect any signal from z > 8 galaxies. The presence of fluctuations 
at optical wavelengths requires a low-redshift signal in addition to high-redshift sources to explain 
combined optical and IR background intensity fluctuations. 
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Figure 7 | Intrahalo light fraction and model intensities (a) shows /ihl, the intrahalo light frac¬ 
tion, as a function of halo mass. The dark and light shaded regions show the 95% and 68% ranges 
of /ihl from anisotropy measurements, and from an analytical predictionpl (blue). Intracluster 
measurements are shown as boxed^, with 1 a errors. The red downward arrows denote the 95% 
confidence upper limit on /ihl estimated for Andromeda (M31) and our Milky Way (MW), fol¬ 
lowing Figure 2 of Ref. (h), d{X I\)/dz from the model, as a function of redshift. We show 
the 68% confidence uncertainties derived from MCMC fitting of the data at 0.606, 0.775, 0.850, 
1.25 and 1.6 um. The total IHL intensity is O.lS/n'n?, 0.24/n 0.28/nO.dS/nn?, and 0.54/n o? 

nW m-2 sr-i for 0.606, 0.775, 0.850, 1.25 and 1.6 /im, respectively. 
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Parameter 

Best fit 

Best fit (no high-x) 

Prior min, max 

loglo(A8<z<13) 

1 iq+<J.27 
-‘-•■‘-^-2.62 

- 

-5,7 

loglo(AiHL) 

q 9q-|-0.14 

'J-^'J-0.12 

q qo-t-0.25 
'J-'J^-0.09 

-6, 10 

a 

-‘-•UU_o.99 

1.00-0.73 

-5,5 

f\ow—z 

0.47 ±0.03 

0.47 ±0.03 

0 .1, 10 

A 1.6 
^DGL 

(3.74+°-30) X 10^ 

(3.72l‘];i) X 10" 

103, 103 

A 1.1&1.25 
^DGL 

(4.3510;57^) X 10^ 

(4.72+iJ:^2) ^ 

103,103 

A 0.850 
^DGL 

(2.83l[5-^^) X 10^ 

(2.771°;32) X 103 

102,10" 

A 0.775 
^DGL 

(2.74ii) X 103 

(2.65l°-33) X 103 

102,10" 

A 0.606 
^DGL 

(1.61«i°) X 103 

(1.43l°:23) X 103 

102, 10" 

pl.6 

'^£,shot 

(7.54 ±0.13) X 10-11 

(7.54 ±0.13) X 10-11 

10-11,10-1° 

pi.25 
'^£,shot 

(7.77l°|i) X 10-11 

(7.77 ±0.14) X 10-11 

10-11, 10-1° 

p0.850 

'^£,shot 

(7.731°;!^) X 10-12 

(8.10 ±0.45) X 10-12 

10-12,10-11 

pO.775 

'^£,shot 

(4.60l°i°) X 10-12 

(4.65 ± 0.30) X 10-12 

10-12, 10-11 

p0.606 

'^£,shot 

(3.271°;21) X 10-12 

(3.39 ±0.15) X 10-12 

10-13, 10-11 


Table 11 Summary of free model parameters. The best-fit values are quoted with la errors. 
logio(A 8 <z<i 3 ) is the high-redshift component used to constrain the SFRD during the reionization 
epoch, which is fit to the 1.25 and 1.60 /rm bands. logio(AiHL) and a are the two parameters 
necessary to describe the IHL component (to wit: in Equation 9). /low-z is the low redshift 

scaling factor which varies the low redshift power spectrum within ala uncertainty. A^ql and 
Qshot respectively the DGL amplitude scaling factor and shot noise at wavelength i. All 
parameter values have units of (nW m“^ sr“^)^. 
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Figure 8 | Probability distributions of fitted model parameters. Here we show the probablity 
density distributions for our fitted model parameters logio(AiHL), a, /low-z, and logio(Ahigh-z) 
eorresponding to the distribution from 8 < 2 ; < 13. The single eurves on the outermost eolumn 
of each row, labeled with a “P”, show the marginalized probability distribution for each parameter 
labeled on the bottom of the figure. Contour regions to the left of these probability distributions 
show how the parameters scale with one another. Each of the shaded regions in the contours 
correspond to the 1, 2 and 3 a uncertainty ranges. 
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Figure 9 | The UV luminosity density and star-formation rate density as measured with inten¬ 
sity fluctuations. Plotted here is the specific UV luminosity density (left axis), with the equivalent 
star formation rate density (SFRD, right axis), as a function of the redshift 2 ;. We show the 1 a and 
2 a error bounds in our redshift bin as the light and dark blue regions. Results from low-redshift 
surveys are shown as blue triangle^, bright green squared, and orange pentagond^. At 2 ; ~ 4 
to 10 the star formation rate density is shown to decrease with increasing redshift as measured by 
previous works, plotted as filled cyan circled, filled red circled^, open red circled filled green 
circled^ll and open blue circled. Gamma ray burst (GRB) studies are shown as gray triangled^, 
squared^ and dark gray circled^. Except for Ref. 1^, other estimates are luminosity function ex¬ 
trapolations and integrations down to Myv = —13. Our measured star formation rate densities are 
consistent with previous works between z ~ 8 to 10, however only extremely bright galaxies are 
directly detected in the aforementioned works with extrapolations down to Muv = —13 involves 
the measured faint-end slope of the luminosity function. For reference we plot the theoretically 
expected relational between UV luminosity density and redshift to reionize the universe and/or to 
maintain reionization using an optical depth to reionization of r = 0.066 ± O.OI 2 P. We take a gas 
clumping factor of G = 3 and show two cases where the escape fraction of galaxies /esc is 6% and 
20 % (see also Ref.l^. 
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Supplementary Figure 1 | Mode-mode coupling correction, a, the mode-mode eoupling matrix 
as generated from our souree mask, b, validity simulation. We generate 90 Gaussian maps 
with a known input power speetrum (blaek filled eireles). For eaeh realization, we apply the mask 
to the simulated map and compute the resulting power spectrum (blue diamonds). Finally, we 
correct the masked power spectrum with our mode-mode coupling matrix to recover the input 
power spectrum (red crosses), which validates our masking correction. 
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Supplementary Figure 2 | Transfer and beam functions. The T(£)’s were generated from a min¬ 
imum of 50 simulations in eaeh band. These simulations ineorporate the effeets of the map-making 
algorithm, tiling pattern, varying exposure depths, residual (temporal) offsets, and eropping effeets, 
speeifie to eaeh filter. The beam transfer funetion, B{i), in eaeh band is just the PSF of eaeh band 
in harmonie spaee. 
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£^ C\^ / 1-K 

e Ci-^s / 27r 

£2 qO. 850 / 27r 

/ i-K 

P (70.606 ! 2t^ 

1.81 X 10=^ 
2.64 X 10® 
4.04 X 10® 
6.40 X 10® 
1.04 X 10-* 
1.71 X 10-* 
2.85 X lO"* 
4.76 X lO"* 
7.99 X 10-* 
1.34 X 10® 
2.26 X 10® 

3.82 X 10® 
6.44 X 10® 
1.09 X 10® 

1.10 ± 1.17 
(8.30 ± 6.60) X 10-1 
(7.78 ±4.09) X 10-1 
(6.33 ± 2.30) X 10-1 
(4.18 ± 0.94) X 10-1 
(4.01 ± 0.63) X 10-1 
(2.89 ± 0.39) X 10-1 
(2.76 ± 0.37) X 10-1 
(2.50 ± 0.43) X 10-1 
(3.19 ± 0.56) X 10-1 
(7.29 ± 0.80) X 10-1 
1.62 ±0.09 

4.59 ±0.13 
(1.49 ± 0.02) X IQi 

1.10 ±0.98 

1.08 ±0.72 
(7.94 ±4.04) X 10-1 
(9.61 ±3.53) X 10-1 
(6.60 ± 1.43) X 10-1 
(5.29 ± 0.77) X 10-1 
(3.71 ± 0.46) X 10-1 
(3.39 ± 0.41) X 10-1 
(2.85 ± 0.46) X 10-1 
(4.03 ± 0.63) X 10-1 
(8.50 ± 0.89) X 10-1 

1.80 ±0.11 

4.96 ±0.16 
(1.55 ± 0.05) X IQi 

(1.86 ± 1.33) X 10-1 
(9.55 ± 5.10) X 10-2 
(8.30 ± 3.25) X 10-2 
(9.61 ±3.23) X 10-2 
(6.70 ± 1.47) X 10-2 
(3.13 ± 0.41) X 10-2 
(2.66 ± 0.28) X 10-2 
(2.90 ± 0.27) X 10-2 
(2.54 ± 0.30) X 10-2 
(3.76 ± 0.45) X 10-2 
(7.76 ± 0.69) X 10-2 
(2.33 ± 0.13) X 10-1 
(5.73 ± 0.48) X 10-1 
2.31 ±0.80 

(2.30 ± 1.69) X 10-1 
(9.02 ± 6.19) X 10-2 
(7.21 ± 3.29) X 10-2 
(9.40 ± 3.34) X 10-2 
(3.98 ± 0.91) X 10-2 
(3.70 ± 0.57) X 10-2 
(3.02 ± 0.39) X 10-2 
(2.30 ± 0.34) X 10-2 
(1.92 ± 0.40) X 10-2 
(2.75 ± 0.55) X 10-2 
(5.16 ± 0.74) X 10-2 
(1.20 ± 0.10) X 10-1 
(3.24 ± 0.22) X 10-1 
1.09 ±0.47 

(1.70 ± 1.30) X 10-1 
(9.01 ± 5.38) X 10-2 
(4.74 ± 2.47) X 10-2 
(3.29 ± 1.17) X 10-2 
(2.71 ±0.57) X 10-2 
(1.86 ± 0.27) X 10-2 
(1.24 ± 0.17) X 10-2 
(9.62 ± 1.74) X 10-^ 
(1.13 ± 0.24) X 10-2 
(1.67 ± 0.34) X 10-2 
(3.27 ± 0.47) X 10-2 
(8.96 ± 0.75) X 10-2 
(2.26 ± 0.13) X 10-1 
(7.08 ± 0.70) X 10-1 


Supplementary Table 1 | Final HST power spectra. Corrected auto-spectra, in units of 

(nW m~^ sr“^)^, for 5 bands. The quoted errors are the 1 a uncertainties. 
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Supplementary Note 1 | Data Availability: The self-calibrated mosaics used for the fluc¬ 
tuation study, including jack-knives with data separated to epochs, and the detected source 
mask, mode-coupling matrix, beam functions, and the transfer functions are available at 
http://herschel.uci.edu/CANDELS 
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